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ABSTRACT 


For  most  applications  in  radar  data  processing,  the  Fourier 
transform  performs  satisfactorily.  However,  other  methods  of 
spectral  analysis  can  offer  some  advantages  when  a data  set  is 
too  short  for  a Fourier  transform  to  resolve  or  detect  important 
spectral  features.  This  report  describes  one  alternative 
technique,  maximum  entropy  spectral  analysis  (MESA),  and  suggests 
possible  radar  applications  including  range-Doppler  sizing  and 
the  coherent  measurement  of  range  rate.  Practical  examples 
demonstrate  an  improvement  in  velocity  resolution  and  cross - 
range  resolution.  Computer  codes  are  listed  that  calculate 
MESA  power  spectra  for  a real  or  complex  discrete  time  series. 
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1. 


INTRODUCTION 


The  interpretation  of  coherent  radar  signals  often  involves 
the  calculation  of  a power  spectrum  or  a periodogram  that 
describes  the  frequency  content  of  data.  The  conventional 
Fourier  approach,  based  on  the  work  of  Wiener  (1950)  and  of 
Blackman  and  Tukey  (1959) , relates  the  autocorrelation  function 
of  a signal  and  its  power  spectrum  through  the  Fourier  transform. 
Cooley  and  Tukey  (1965)  popularized  the  Fourier  approach  with 
the  computationally  efficient  fast  Fourier  transform  (FFT) , 
which  has  dominated  the  analysis  of  radar  data  until  the  present 
time . 

Notwithstanding  the  speed  and  mathematical  elegance  of  the 
Fourier  transform  with  its  entourage  of  weighting  functions  and 
tapering  schemes,  it  is  troubled  by  several  unavoidable  limita- 
tions that  become  serious  as  the  temporal  length  of  the  data 
set  shortens.  However,  alternative  methods  of  spectral  analysis 
do  exist  that  have  not  as  yet  found  wide  application  in  radar. 
Among  these  techniques  is  "maximum  entropy  spectral  analysis" 
(MESA),  an  outgrowth  of  the  predictive  deconvolution  filtering 
techniques  developed  by  geophysicists  for  oil  exploration.  The 
primary  purpose  of  this  report  is  to  introduce  MESA  to  radar  as 
a means  of  improving  the  velocity  resolution  and  cross-range 
resolution  currently  limited  by  conventional  Fourier  concepts. 


Section  2 briefly  identifies  some  of  the  problems  with  the 
Fourier  approach  to  spectral  analysis.  Section  3 describes  the 
MESA  procedure  for  which  more  detailed  derivations  are  given  in 
the  appendices.  Comparing  conventional  and  MESA  spectra, 
Section  4 treats  several  radar  applications  to  include  the 
measurement  of  range  rate  and  the  range-Doppler  sizing  of  hard 
bodies.  Many  other  applications  are  feasible.  In  principle, 
MESA  can  replace  a Fourier  transform  wherever  the  latter  may 
occur,  whether  a transformation  is  to  be  made  from  the  time 
domain  to  the  frequency  domain,  or  from  the  frequency  domain  to 
the  time  domain. 
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2.  CONVENTIONAL  FOURIER  TECHNIQUES 

Conventional  methods  of  spectral  analysis  make  unrealistic 
assumptions  about  the  nature  of  data  outside  of  the  observation 
interval.  The  periodogram  approach  (Jones,  1965)  projects  the 
data  set  as  an  infinite  periodic  repetition  of  itself,  while 
the  familiar  autocorrelation  approach  (Blackman  and  Tukey, 

1959)  assumes  a zero  extension.  These  schemes  attempt  to 
lessen  the  effects  of  the  finite  length  of  a data  set  upon 
Fourier  transformation. 

One  may  think  of  a data  set  truncated  in  the  time  domain 
as  the  product  of  an  infinite  data  set  and  a window  function 
that  has  non- zero  amplitude  only  within  the  observation  interval. 
The  measured  power  spectrum  of  this  product  is  the  convolution 
in  the  frequency  domain  of  the  true  spectrum  of  the  data  set 
(which  one  would  like  to  have  measured)  and  the  spectrum  of  the 
window  function  itself.  The  measured  spectrum  is  unavoidably  a 
distortion  of  the  true  spectrum;  the  convolution  allows  the 
calculation  of  the  spectrum  at  one  frequency  to  be  contaminated 
by  energy  at  all  other  frequencies  (loosely  termed  "leakage 
through  the  sidelobes"  of  the  window  function).  The  measured 
power  spectrum  can  even  take  on  physically  meaningless  negative 
values  if  negative  lobes  in  the  window  spectrum  are  being 
convolved  with  strong  frequency  components  in  the  true  spectrum. 


For  example,  a rectangular  window  (the  default  window)  has 
a spectrum  of  the  form 

W(f)  - sin(nfT)/itfT 

where  T is  the  length  of  the  data  set.  The  first  zero  occurs 
at  f ■ 1/T  and,  by  rule  of  thumb,  frequency  components  in  the 
true  spectrum  spaced  more  closely  than  1/T  are  not  easily 
resolved  in  the  measured  spectrum  owing  to  the  smearing  caused 
by  convolution.  The  frequency  resolution  of  the  measured 
spectrum  is  therefore  6f  •»  1/T. 

In  radar  applications,  frequency  shift  and  velocity  are 
related  by  the  equation 

v ■ Xf/2 

so  that  velocity  resolution  becomes  the  familiar  result 
6v  - X6f/2  - X/2T. 

Here  X is  the  radar  wavelength.  We  observe  that  this  limitation 
on  the  resolution  of  velocity  is  a result  of  using  conventional 
Fourier  spectral  analysis  and  does  not  necessarily  apply  if 
another  method  of  spectral  analysis  is  used. 

The  rectangular  window  is  sometimes  replaced  by  a tapering 
function  (Hamming,  Bartlett,  and  Taylor  windows  are  a few 
examples),  designed  to  reduce  the  effect  of  the  window  sidelobes 
when  the  window  spectrum  is  convolved  with  the  true  spectrum. 
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The  choice  of  a window  always  involves  a compromise  between 
frequency  resolution  and  the  extent  to  which  window  sidelobes 
allow  the  calculation  of  the  measured  spectrum  at  one  frequency 
to  be  contaminated  by  components  at  all  other  frequencies. 

The  Fourier  transform  loses  resolution  when  the  data  set 
is  short  compared  to  the  periods  of  the  spectral  components  in 
the  data.  Toman  (196S)  and  Jackson  (1967)  have  shown,  for 
example,  that  most  of  the  energy  of  a sinusoid  will  appear  near 
zero  frequency  in  a spectrum  calculated  with  a Fourier  transform 
if  the  length  of  the  data  set  is  less  than  581  of  the  period  of 
that  sinusoid.  Interference  effects  caused  by  the  window 
sidelobes  can  also  produce  spectral  shifts  as  large  as  161. 

Basically,  the  problems  from  which  conventional  techniques 
suffer  are  caused  by  the  finite  length  of  the  data  set.  While 
window  theory  is  elegant,  it  treats  only  the  symptoms  of  trunca- 
tion and,  in  doing  so,  corrupts  perfectly  good  data  with  weight- 
ing functions.  Abies  (1974)  cites  what  he  calls  the  "First 
Principle  of  Data  Reduction": 

The  result  of  any  transformation  imposed 
on  experimental  data  shall  incorporate 
and  be  consistent  with  all  relevant 
data,  and  be  maximally  non-committal 
with  regard  to  unavailable  data. 

Conventional  Fourier  techniques  clearly  violate  this 
principle:  data  not  in  evidence  are  unrealistically  assumed 
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and  true  data  are  changed  by  weighting  functions.  Conventional 
techniques  are  therefore  philosophically  objectionable  and 
become  less  reliable  when  the  data  set  becomes  short  - and  it 
is  precisely  the  short  data  set  with  which  radar  is  usually 
concerned. 
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3. 


MAXIMUM  ENTROPY  SPECTRAL  ANALYSIS  (MESA) 


3.1  Deconvolution  Filtering 

The  maximum  entropy  approach  to  spectral  analysis  is  a 
variation  of  the  deconvolution  filtering  techniques  developed 
by  geophysicists  for  processing  seismic  signals.  A deconvolu- 
tion filter  whitens  the  spectrum  of  the  signal  on  which  it 
operates;  that  is,  when  convolved  with  the  original  signal,  it 
outputs  a new  signal  with  a constant  (white)  spectrum.  Mathemat- 
ically expressed,  the  convolution  of  a discrete-time  series 
x(n)  with  a digital  filter  with  coefficients  a(n)  is 

M 

f(n)  * £ x(n)  a(n  - k) . (3-1) 

k»o 

The  spectrum  of  f(n),  which  is  the  product  in  frequency  space 
of  the  spectrum  of  x(n)  and  the  transfer  function  (or  impulse 
response  function)  A(f)  of  the  filter,  is  a constant: 

F (f)  - constant  - X(f)  A(f).  (3-2) 

The  power  spectrum  of  x(n)  is  then 

P(f)  - |X(f)|2  - K/ 1 A(f) | 2 (3-3) 

where  K is  a constant.  Simply  stated,  deconvolution  filtering 
involves  finding  the  digital  filter  that  changes  the  input 
signal  into  an  output  signal  with  a constant  spectrum. 
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This  approach  to  spectral  analysis  is  also  known  as  the 
Markov  spectrum  or  the  autoregressive  spectrum  described  by 
statisticians  (see,  for  example,  Parzen,  1969).  Burg  (1967; 

1968;  1975)  realized  that  this  approach  yields  the  spectrum 
having  the  "maximum  entropy"  (explained  below)  of  all  possible 
spectra  that  are  consistent  with  the  measured  autocorrelation 
function  of  x(n).  Burg  also  devised  methods  of  efficiently 
calculating  the  filter  coefficients  from  which  A(f)  can  be 
determined. 

One  advantage  of  deconvolution  filtering  is  immediately 
obvious:  finding  X(f)  does  not  involve  a convolution  in  frequen- 

cy space  with  a cumbersome  window  spectrum  that  (unavoidably) 
destroys  spectral  resolution.  The  convolution  has  already 
taken  place  in  the  time  domain  between  the  input  signal  and  the 
digital  filter,  Therefore,  there  arc  no  window  sidelobes  or 
serious  end  effects  as  occur  upon  conventional  Fourier  transfor- 
mation. The  truncation  of  the  data  set  is  important  only  to 
the  extent  that  enough  data  must  be  available  to  allow  the 
construction  of  an  efficient  whitening  filter  that  can  reduce 
the  data  to  a random  series. 

Other  comparisons  between  the  characteristics  of  Fourier 
and  MESA  spectra  will  be  drawn  later  and  are  summarized  in  the 
appendices.  However,  a discussion  of  the  philosophy  of  maximum 
entropy  is  first  in  order. 
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3.2  The  Maximum  Entropy  Philosophy 

Choosing  the  best  spectral-domain  representation  of  a 
truncated  discrete  time  series,  for  which  only  an  imperfectly 
determined  autocorrelation  function  can  be  calculated,  is  a 
major  problem  in  signal  analysis.  Among  the  countless  spectra 
that  may  be  consistent  with  a given  autocorrelation  function, 
only  one  spectrum  can  be  optimal.  A set  of  rules  governing 
that  choice  must  be  established. 

Jaynes  (1957)  introduced  a method  of  statistical  inference 
called  the  "maximum  entropy  estimate".  He  showed  that  informa- 
tion theory  (Shannon  and  Weaver,  1949)  provides  a criterion  for 
selecting  the  best  statistical  description  of  a process  when 
only  a partial  knowledge  of  that  process  is  available.  The 
optimal  choice  is  the  one  which  is  maximally  non-committal  with 
regard  to  any  missing  data,  and  which  is  simultaneously  con- 
strained to  be  consistent  with  all  available  data.  The  result 
is  the  best  estimate  that  could  have  been  made  on  the  basis  of 
the  data  at  hand.  (Actually,  Abies 's  first  principle  is  a 
restatement  of  Jaynes's  conclusion.) 

The  term  "entropy"  is  used  in  an  informational  sense,  such 
that  a measure  of  the  entropy  of  a process  is  a measure  of  our 
ease  in  coping  with  data  that  we  do  not  have.  We  must  make  no 
assumptions  and  impose  no  constraints  that  cannot  be  justified 
directly  with  the  data  that  we  do  have. 


What  is  needed  is  an  expression  for  the  entropy  of  a 
process  in  terms  of  its  spectrum.  That  expression  is  then 
maximized  subject  to  the  constraints  imposed  on  the  spectrum  by 
the  available  data. 

Shannon  and  Weaver  (1949)  have  shown  that  the  appropriate 
expression  for  the  entropy  of  a process  in  terms  of  its  power 
spectrum  P(f)  is 

oo 

E - / log  P(f)  df.  (3-4) 

— OO 

The  values  of  the  autocorrelation  function  calculated  from  the 
available  data  comprise  the  constraints  on  P(f).  Ideally, 
every  value  <j>  (m)  of  the  autocorrelation  function  calculated  at 
lag  m is  related  to  P(f)  by 

OO 

<Km)  * / P(f)  exp  [- i2irmf At ] df  (3-5) 

— oo 

where  At  is  the  time  spacing  of  the  data. 

Maximizing  Eq.  (3-4)  subject  to  the  constraints  of  Eq.  (3- 
5)  (one  for  each  known  value  of  Km))  becomes  a problem  in  the 
calculus  of  variations.  In  the  case  where  the  Km)  are  equally 
spaced  and  centered  at  zero  lag,  P(f)  may  be  determined  with 
Lagrange  multipliers  (Edwards  and  Fitelson,  1973). 
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However,  a more  direct  approach,  which  is  equivalent  to 
the  formal  use  of  Lagrange  multipliers  but  which  allows  a more 
tangible  appreciation  of  the  MESA  procedure,  is  the  determina- 
tion of  the  deconvolution  digital  filter  that  transforms  a 
given  time  series  into  a random  series  with  a white  spectrum. 
The  next  section  and  Appendix  I outline  this  calculation  from 
which  A(f)  in  Eqs.  (3-2)  and  (3-3)  is  found. 

3.3  Prediction-Error  Filtering 

The  whitening  of  a discrete- time  series  can  be  done  with  a 
prediction-error  filter  (Peacock  and  Treitle,  1969;  Makhoul , 
1975).  A complex  measurement  x(n)  is  approximated  by  the 
weighted  average  of  the  preceding  M terms,  or  by  the  weighted 
average  of  the  successive  M terms.  The  former  is  a forward 
prediction,  and  the  latter  is  a backward  prediction: 

M 

xf(n)  = £ aM  . x(n-k)  (3-6) 

f k-1  ’ 

M * 

xfe(n)  = EaM,kx(n+lc^  (3-7) 

k — l 

where  the  asterisk  denotes  complex  conjugation,  and  a^  ^ is  the 
kth  filter  coefficient  out  of  a total  of  M coefficients.  The 
prediction  filter  is  a linear  combination  of  M data  points  that 
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predicts  the  data  point  immediately  preceding  (or  following) 
the  M points.  The  errors  in  the  forward  or  backward  prediction 
of  a given  data  point  x(n)  are 

ef(n)  = x(n)  - xf(n)  (3-8) 

eb(n)  = x (n)  - xfe(n)  (3-9) 

and,  for  the  linear  prediction  to  be  optimal,  the  errors  should 
be  minimized  simultaneously  in  a least  squares  sense.  Therefore, 
we  minimize  the  total  error  power 

PM  ■ £ [efCn)2  * ebfn)2]  C3’™) 

with  respect  to  the  M coefficients.  PM  is  the  power  of  the 
output  error  series  left  behind  after  the  x(n)  have  been  predic- 
ted. The  spectrum  of  the  output  series  is  a constant  because, 
if  the  spectrum  contained  any  recognizable  frequency  components, 
then  we  could  use  the  knowledge  of  those  components  to  improve 
the  prediction  by  editing  the  filter  coefficients  until  the 
spectrum  does  become  white. 

Therefore,  the  prediction  filter  can  be  used  to  transform 
the  data  x(n)  into  another  series  (the  eAror  series)  that  has  a 
constant  spectrum:  each  point  x(n)  is  replaced  by  the  error  in 

its  prediction.  The  predict ion -error  filter  is  represented  by 


12 


Eqs.  (3-8)  and  (3-9)  which  use  the  prediction  filter.  The 
output  of  the  former  is  the  error  in  the  prediction  of  a known 
measurement;  the  output  of  the  latter  is  simply  the  prediction 
itself. 

Van  Den  Bos  (1971)  has  shown  that  the  prediction-error  (P- 
E)  filter  is  equivalent  to  the  least- squares  fitting  of  a 
discrete- time  all-pole  model  to  the  data.  Anderson  (1974)  ex- 
tended this  work  in  real  form,  and  a complex  formulation  of 
Anderson's  algorithms,  which  calculate  the  P-E  filter  coeffi- 
cients, is  derived  in  Appendix  I. 

The  impulse  response  function  of  the  P-E  filter  is  its  z- 
transform 


A(f)  « 1 + £ aM>k  z'K  (3-11) 

k=  1 * 

z = exp  [i2irfAt].  (3-12) 

There  is  one  coefficient  a^  k for  each  of  the  M constraints  put 
on  the  power  spectrum  during  entropy  maximization. 

The  choice  of  M is  somewhat  arbitrary,  although  Akaike 
(1969a,  b;  1970)  has  argued  that  the  length  M of  the  filter 
should  be  chosen  so  as  to  minimize  the  error  power  P^.  An  M 
should  be  chosen  so  that  increasing  the  filter  length  to  M + 1 
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longer  significantly  reduces  the  error  in  prediction;  that 
is,  PM+1  is  not  much  smaller  than  P^.  Of  course,  M cannot 
exceed  the  number  of  data  points.  (See  the  review  article  by 
Ulrych  and  Bishop,  1975.)  Using  too  small  an  M results  in  a 
highly  smoothed  spectrum,  obviating  the  improved  resolution 
capabilities  of  MESA;  using  too  long  a filter  allows  noise  to 
introduce  spurious  detail  into  the  spectrum.  A reasonable 
initial  value  of  M is  one  fifth  the  number  of  the  data  points 
of  the  input  sequence. 

3.4  Power  Spectrum  and  Power  Spectral  Density 

Both  Lacoss  (1971)  and  Burg  (1975)  have  pointed  out  that 
the  MESA  power  spectral  density  function  is  a better  measure  of 
the  amount  of  power  in  a small  bandwidth  than  is  the  power 
spectrum.  The  two  are  related  by 

f*6f/2 

PSD (f ) = / P(f)  df/6f  (3-13) 

f-<5f/2 

where  PSD(f)  is  the  power  spectral  density,  P(f)  is  the  power 
spectrum  (Eq.  (3-3)),  and  <5f  is  the  small  bandwidth  over  which 
P(f)  is  integrated.  Since  P(f)  is  calculable  at  any  frequency, 
A(f)  being  an  analytic  and  continuous  complex  polynomial, 

Eq.  (3-13)  is  easily  implemented  using  numerical  integration 
(see  appendices) . 


3.5  MESA  and  Non-Linearity 


One  attractive  property  of  conventional  Fourier  techniques 
is  linearity:  the  power  spectral  density  of  the  sum  of  several 

signals  is  the  sum  of  the  power  spectral  densities  of  the 
individual  signals.  Moreover,  the  square  root  of  the  Fourier 
PSD(f)  gives  a reliable  amplitude  spectrum  of  the  input  signal. 

Unfortunately,  MESA  is  ultimately  not  a linear  technique. 
Indeed,  Eq.  (3-3)  is  the  inverse  of  the  square  of  a complex 
polynomial . Therefore,  a comparison  of  relative  power  at 
different  frequencies  can  be  misleading.  Parseval’s  theorem 
notwithstanding,  the  integral  of  the  MESA  PSD(f)  may  not  equal 
the  total  power  of  the  input  signal. 

The  value  of  MESA  lies  in  frequency  detection  and  frequency 
resolution,  as  discussed  later  in  Section  3.7.  It  is  not  as 
useful  for  determining  the  relative  strengths  of  different 
frequency  components.  MESA  can  be  used  in  conjunction  with 
Fourier  techniques,  for  example,  by  identifying  the  important 
frequency  components  with  MESA  and  then  performing  a DFT  (dis- 
crete Fourier  transform)  only  at  those  frequencies  to  estimate 
their  strengths  and  phases.  Section  3.6  offers  another  alterna- 
tive in  which  a Fourier  transform  can  be  taken  for  a data  set 
extended  beyond  the  observation  period  by  the  linear  prediction 
filter. 
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It  is  worthwhile  at  this  point  to  mention  the  non-linear 
effects  of  noise  on  the  MESA  spectrum.  We  can  rewrite  Eq.  (3- 
3)  as 

P(f)  - K/  | z"1-  ZjJ  2 | z ~ 1 - z 2 1 2 ...  K1-^2  (3-14) 

z * exp  (iO)  » exp  (i2nfAt) 

where  we  have  factored  A(f)  in  the  denominator.  If  we  plot  the 
poles  z | (i-1  to  M)  of  P(f)  in  the  complex  plane,  as  demonstrat- 
ed in  Fig.  3-1  for  M*5,  they  will  all  lie  within  the  unit 
circle.  (This  is  because  the  set  of  prediction  filter  coeffi- 
cients is  "minimum  phase",  as  explained,  for  example,  in  Robinson 
(1967)  or  Claerbout  (1976).)  As  A(f)  is  evaluated  on  the  unit 
circle,  the  angle  0=27rfAt  varies  from  -tt  to  tt.  Therefore,  f 
varies  between  — l/2At  and  l/2At,  that  is,  between  minus  and 
plus  the  Nyquist  frequency.  For  a given  frequency  at  point  F, 
each  factor  in  the  denominator  of  Eq.  (3-14)  is  the  square  of 
the  distance  between  point  F and  one  of  the  poles.  As  f changes 
and  point  F passes  by  a pole  near  the  unit  circle,  P(f)  will 
exhibit  a local  maximum. 

The  frequency  components  in  the  input  data  set  will  corre- 
spond to  those  poles  of  P(f)  closest  to  the  unit  circle. 

Clearly,  there  can  be  no  more  maxima  detected  in  P(f)  than 
there  are  poles  of  P(f)  or  coefficients  of  the  filter.  In  the 
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case  where  there  are  actually  fewer  frequency  components  in  the 
data  than  there  are  poles  available  for  placement , the  extra 
poles  can  be  eliminated  by  reducing  the  number  of  filter  coeffi- 
cients, as  M is  the  degree  of  the  complex  polynomial  A ( f ) that 
is  factored  to  locate  he  M poles. 

Noise  can  affect  the  location  of  the  poles  in  the  complex 
plane.  Since  Fq.  (3- 14 ) involves  the  inverse  squares  of  the 
distances  between  a point  on  the  unit  circle  and  the  poles,  a 
small  change  in  the  position  of  a pole  that  is  already  near  the 
unit  circle  can  cause  a large  change  in  the  magnitude  of  P(f) 
near  that  pole.  It  is  for  this  reason  that  the  amplitude  of 
P(f)  may  not  accurately  reflect  the  true  power  spectrum.  Also, 
since  the  PSD(f)  is  the  line  integral  of  P(f)  over  a small  arc 
on  the  unit  circle  (a  small  frequency  interval  <5f),  the  ampli- 
tude of  the  power  spectral  density  similarly  may  not  be  accurate. 

Burg  (197b)  has  recognized  that  both  the  peak  value  and 
width  of  an  apparent  spectral  line  in  P(f)  strongly  depend  on 
the  background  noise.  However,  he  maintains  th?t  the  product 
of  the  peak  value  and  line  width,  which  is  proportional  to  the 
total  power  in  the  spectral  line,  will  be  estimated  accurately 
if  the  sampling  of  the  spectrum  is  fine  enough  lo  trace  out  the 
shape  of  the  spectrum.  In  general,  the  unequivocal  detection 
of  a specttal  component  is  much  more  reliable  than  an  estimate 


17 


of  its  strength  in  a MESA  spectrum.  Therefore,  locating  the 
poles  of  P(f)  may  be  sufficient  for  frequency  detection. 

Akaike  (1969),  Baggeroer  (1976),  and  others  cited  therein 
have  investigated  some  confidence  intervals  and  the  statistical 
variability  of  the  MESA  technique. 

The  MESA  user  must  decide  if  the  non-linearity  poses  a 
problem  for  his  specific  application.  If  the  detection  of 
frequency  components  in  a data  set  is  the  primary  goal,  MESA 
will  indeed  be  useful.  However,  the  procedure  may  not  be 
totally  adequate  if  precise  amplitude  and  phase  information  is 
also  required.  In  that  case,  a DFT  can  be  computed  at  the 
frequencies  already  identified  by  the  MESA  procedure. 

3.6  Linear  Prediction  and  the  Fourier  Transform 

The  MESA  procedure  is  linear  up  to  the  point  at  which  the 
prediction  filter  is  calculated.  The  filter  itself  is  a linear 
operator  that,  as  we  have  seen,  uses  a weighted  average  of  M 
data  points  to  estimate  the  adjacent  data  points. 

We  propose  the  following  procedure  to  take  advantage  of 
the  linearity  of  both  the  prediction  filter  and  the  Fourier 
transform:  A prediction  filter  is  calculated  from  a data  set 

of  points.  The  filter  is  then  used  to  extend  (or  predict) 
the  original  data  set  to  a total  of  N2  points,  outside  of  the 
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observation  interval  in  both  forward  and  backward  directions. 

Then  u conventional  Fourier  spectrum  is  calculated  with  the 
larger  data  set  of  N2  points. 

The  data  points  that  are  predicted  by  the  filter  will  have 
essentially  the  same  spectral  composition  as  the  original  data, 
as  a knowledge  of  the  spectrum  of  the  original  data  is  contained 
in  the  filter  coefficients. 

Since  additional  signal  is  being  created  and  subsequently 
transformed,  the  normalization  of  the  power  spectrum  of  the 
extended  data  set  should  be  done  with  respect  to  the  total 
number  of  samples  N2«  Any  weighting  window  should  be  applied 
across  the  N2  samples,  after  the  extension  has  been  done. 

The  deterministic  spectral  components  having  been  reinforced 
in  the  process  of  prediction,  the  apparent  signal- to-noise 
ratio  is  improved;  the  contribution  of  noise  cannot  be  predicted 
and  therefore  does  not  appear  in  the  extension.  Of  course,  the 
quality  of  the  prediction  will  depend  on  the  extricable  determin- 
ism and  the  SNR  of  the  original  data.  If  the  original  data  are 
pure  noise,  the  poles  of  P(f)  will  be  near  zero  (clustered  at 
the  origin  of  the  complex  plane)  and  the  magnitudes  of  the 
predicted  points  will  be  very  small  or  zero.  Otherwise,  the 
sinusoidal  components  of  the  original  data  dominate  the  extension. 


Linear  preoiction  is  less  likely  to  improve  the  regions  of 
a power  spectrum  for  which  a band  or  continuum  of  frequencies 
contains  energy.  In  this  case,  discrete  poles  will  position 
themselves  in  the  complex  plane  in  the  vicinity  of  the  arc  on 
the  unit  circle  corresponding  to  that  band  of  frequencies. 

Some  frequencies  in  the  band  may  be  enhanced  slightly  more  than 
others  because  the  poles  are  points,  each  of  which  is  closest 
to  a single  point  (frequency)  on  the  unit  circle.  The  uneven 
effect  will  be  small,  however,  and  the  poles  will  move  farther 
away  from  the  arc  on  the  unit  circle  as  the  band  of  frequencies 
broadens.  Indeed,  if  there  is  equal  energy  at  all  frequencies 
(noise),  the  poles  cluster  near  the  origin,  which  is  the  only 
point  equidistant  from  every  point  on  the  unit  circle. 

Predicting  the  data  beyond  the  observation  interval  is 
more  palatable  than  assuming  such  data  are  zeroes.  The  predic- 
tion is  as  self-consistent  with  the  original  data  as  is  possible 
with  the  information  at  hand.  Conventional  zeroes  would  be  no 
more  consistent  with  the  original  data  as  would  be  any  other 
constant  randomly  chosen. 

The  appendices  contain  a FORTRAN  program  "LNPRED"  that 
linearly  extends  a complex  data  set  from  to  N2  points.  A 
user- supplied  Fourier  transform  can  then  be  performed  on  the  N2 
points  output  from  this  prediction  subroutine. 
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3.7  Improved  Frequency  Resolution  with  MESA 

Because  a data  set  could,  in  principle,  be  extended  to 
infinity,  spectra  based  on  linear  prediction  have  a better 
chance  both  of  separating  closely  spaced  frequency  components 
and  of  locating  them  more  precisely  in  frequency  space.  In 
this  section,  we  demonstrate  that  MESA  can  more  accurately 
estimate  the  frequency  of  a single  noisy  sinusoid,  as  well  as 
detect  two  closely  spaced  sinusoids,  than  can  conventional 
techniques . 


First  consider  a noisy  sinusoid  sampled  in  time: 

» 

x(nAt)  ■ A sin(2irfonAt  * $)  + N(nAt)  (3-15) 

where  fQ  is  the  basic  frequency,  $ is  a phase  constant.  At  is 

the  sample  spacing,  A is  the  sinusoid  amplitude,  n is  a sample 

counter,  and  N is  Gaussian  noise.  We  want  to  calculate  the  [ 

error  in  measuring  fQ , as  a function  of  a given  signal-to-noise 

power  ratio  (A/N)  and  of  a fixed  number  of  data  samples,  when  * 

i: 

either  MESA  or  conventional  techniques  are  used. 

••a 

<1 

We  assemble  an  ensemble  of  100  data  sets  {x(nAt)}  by  | 

making  100  random  draws  on  <t>  between  0 and  2w.  Each  data  set 
has  the  same  number  of  samples  and  the  same  signal  and  noise  p 

amplitudes.  A spectral  transform  is  done  on  each  of  the  100 
data  sets  in  an  attempt  to  extract  fQ.  Then  the  standard 
deviation  of  the  100-member  set  ^measure(j  — fQi  is  calculated, 
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and  cr/ f 0 becomes  a measure  of  the  precision  with  which  we  can 
find  f for  a given  SNR  and  number  of  samples.  We  can  con- 
struct a family  of  curves  by  holding  the  SNR  constant  and 
varying  the  number  of  samples,  or  vice  versa. 

As  long  as  f < 0.5/At,  the  sinusoid  is  properly  sampled 

and  the  actual  values  of  f and  At  are  irrelevant. 

o 

Figures  3-2 (a-b)  compare  the  expected  error  in  measuring 
f for  a data  set  that  is  processed  with  either  MESA  or  an  FFT. 
We  have  taken  f - 0.25/At  (four  samples  per  period). 

As  expected,  either  an  improvement  in  SNR  or  an  increase 

in  the  number  of  samples  reduces  the  error.  However,  for  a 

given  SNR  and  length  of  data  set,  the  error  in  the  MESA  estimate 

of  fQ  is  consistently  smaller  than  the  error  in  the  FFT  estimate 

of  f . 
o 

This  reduction  in  expected  error  demonstrates  that  conven- 
tional Fourier  techniques  do  lose  some  information  owing  to  the 
problems  discussed  in  Section  2.  Although  we  do  not  expect  the 
Fourier  transiorm  of  a linearly  predicted  data  set  to  render  an 
error  as  small  as  MESA,  it  should  show  improvement  over  the 
conventional  Fourier  transform  that  uses  an  extension  of  zeroes. 

As  an  example  of  the  ability  of  MESA  and  Fourier  tech- 
niques to  resolve  closely  spaced  frequencies,  consider  a two- 
component  signal  sampled  in  time 


2 

x(nAt)  = £ A^fcos  (2ir£mnAt)  + i s in (2iTfmnAt)  ] + N(nAt) 

m=l 

(3-16) 

The  shaded  part  of  Fig.  3-3  displays  25  samples  of  the 
real  part  of  Eq.  (3-16),  for  which  A^=A2=  2 units,  At  = 0.01  sec, 
f-^  = -15.3  Hz,  and  f2  = -13.3  Hz.  Enough  white  noise  is  added 
to  give  a SNR  of  10  dB. 

Since  only  0.25  sec  of  data  is  transformed,  conventional 
Fourier  techniques  (zero  extension)  will  be  unable  to  resolve 
frequencies  spaced  more  closely  than  1/0.25  = 4 Hz.  In  Fig.  3- 
4 the  conventional  power  spectral  density  shows  only  one  peak, 
as  expected.  The  MESA  PSD  of  Fig.  3-5,  however,  detects  the 
presence  of  both  frequencies  in  the  25  samples;  the  unequal 
power  amplitudes  are  the  result  of  non-linearity. 

If  we  extend  the  25  samples  to  100  samples  with  a 5-point 
prediction  filter  (the  extension  is  the  unshaded  region  of 
Fig.  3-3),  then  the  conventional  Fourier  PSD  does  resolve  both 
frequencies,  as  is  shown  in  Fig.  3-6.  The  power  amplitudes  in 
Fig.  3-6  are  almost  equal  and  yield  a more  reliable  estimate  of 
the  relative  strengths  of  the  two  frequency  components. 

Thus,  we  find  that  both  the  MESA  spectrum  of  the  original 
data  set  and  the  Fourier  spectrum  of  the  linearly-predicted 
data  set  can  resolve  closely  spaced  frequencies  with  more 
success  than  can  the  traditional  Fourier  spectrum. 


23 


<*> 'OkMuiU  AU*.'.  i 


*-»  liim  iMawwi 


mmmirn 


r~" 

% 


3.  8 Additional  Use  of  the  z-Transform 

The  calculation  of  a power  spectrum  or  a power  spectral 
density  is  really  unnecessary  if  only  frequency  information  is 
required.  The  z-tTansform  A(f)  of  the  prediction  filter  can  be 
factored  to  give 

A(f)  = (z'1-z1)(z"1-z2)  ...  (z_1-zM)  (3-17) 

as  in  Eq.  (3-14),  where  z ^(£  = 1 to  M)  are  the  zeroes  of  A(f) 
and  z = exp(i0)  = exp  ( i 2 tt f At) . The  magnitudes  of  the  zeroes 
| | are  always  between  zero  and  unity,  and  the  zeroes  with  the 
largest  magnitudes  are  likely  to  correspond  to  the  frequency 
components  in  the  data.  If  we  choose  some  minimum  magnitude 
which  z ^ must  have  to  be  considered  as  a candidate  frequency 
component  ( | j > 0.8,  for  example),  then  we  can  calculate  the 
corresponding  frequency  since 

0 = 2-rrfAt  = arctan[Im(z|  )/Re  (z^  ) ] . (3-18) 

For  radar  applications,  range  rate  R = Xf/2  so  that 

R = Vamb  arctan[Im(z£)/Re(Zjji)]/2iT  (3-19) 

where  Vamb  = A/2At  is  the  ambiguous  velocity.  Equation  (3-19) 
is  easily  interpreted  geometrically  because  one  trip  around  the 
unit  circle  is  one  foldover  in  velocity.  A target's  range  rate 
can  be  estimated  accurately  in  this  way  without  the  use  of  a 
Fourier  transform. 
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3.9  Additional  Comments 


This  section  has  presented  the  essence  of  the  linear 
prediction  - maximum  entropy  approach  to  spectral  analysis, 
avoiding  many  of  the  mathematical  complexities  that  are  treated 
extensively  in  the  literature  and  briefly  here  in  the  appendices. 
The  reader  interested  in  more  detail  is  referred  to  the  biblio- 
graphy in  the  review  article  by  Ulrych  and  Bishop  (1975)  and  to 
Burg's  (1975)  dissertation  which  explores  many  topics  not 
discussed  here. 
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Fig. 3-1.  Zeroes  of  the  z-transform  of  a five-point 
prediction  error  filter  are  plotted  in  the  complex  plane 
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Fig. 3-3.  A two- frequency  signal  (-15.3  Hz  and  -13.3  Hz) 
is  sampled  25  times,  every  0.01  second  (shaded  region), 
and  linearly  predicted  to  a total  of  100  samples. 
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Fig, 3-4.  Conventional  Fourier  PSD  of  shaded  region  in  Fig. 3-3 
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Fig. 3-5.  MESA  PSD  of  shaded  region  in  Fig. 3 -3. 

I: 

■ 


31 


4. 


MESA  APPLICATIONS  FOR  RADAR 


4.1  Radar  Data  Preparation 

A coherent  radar  records  the  amplitude  and  phase  of  the 
energy  coming  from  each  range  cell  along  the  radar  line  of 
sight  (RLOS) . For  a given  range  cell  or  range  gate,  we  form  a 
complex  sample  at  time  t 

x(t)  ■ A(t)  cos0(t)  + iA(t)  sinO(t)  (4*1) 

where  A(t)  is  the  radar  cross  section  in  volts  (or  in  meters  or 
an  equivalent  linear  unit)  and  6(t)  is  the  phase.  Since  the 
data  are  recorded  pulse  by  pulse,  x(t)  is  a discrete- time 
series. 

If  the  object  is  moving,  then  0(t)  changes  in  time  since 
phase  is  governed  by 

0(t)  • 4irR(t)/X,  (4-2) 

where  R(t)  is  the  one-way  distance  between  the  object  and  the 
radar,  and  X is  the  radar  wavelength.  Pulses  must  occur  often 
enough  that  0(t)  does  not  change  more  than  2ir  between  pulses; 
otherwise,  the  phase  is  aliased  and  becomes  ambiguous. 

The  amplitude  A(t)  may  change  in  time  if  the  orientation 
of  the  radiation  pattern  of  the  object  changes  with  respect  to 
the  RLOS,  owing  to  the  overall  body  velocity  or  to  body  motion 
about  the  center  of  mass  (spin,  tumble,  precession,  etc.). 
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Hence,  x(t)  contains  information  about  the  scintillation 
of  reflected  energy  and  about  the  velocities  of  objects  in  the 
range  cell.  Indeed,  A(t)  and  0(t)  represent  the  resultant 
vector  sum  of  the  individual  returns  from  all  unresolved  scatter- 
ers  in  the  same  range  cell,  each  of  which  contributes  a frequency 
component  to  A(t),  to  0(t),  or  to  both. 

A non-coherent  radar  does  not  record  phase,  for  which  case 
we  take  0(t)  ■ 0 in  Eq.  (4-1)  and  x(t)  reduces  to  the  real 
series 

x (t ) « A(t)  + iO  = A(t ) . (4-3) 

In  this  section,  we  present  examples  of  the  use  of  Fourier 
and  MESA  algorithms  to  measure  body  velocity  (range  rate)  and 
to  perform  range- Doppler  sizing  of  an  object  with  motion  about 
its  center  of  mass. 

4.2  Measurement  of  Range  Rate 

The  projection  of  the  velocity  vector  of  an  object  along 
the  RLOS  is  its  range  rate,  which  can  be  estimated  from  a 
coherent  data  set  by  the  spectral  analysis  of  x(t).  For  the 
moment  we  assume  A(t)  is  a constant  and  temporal  changes  in 
x(t)  can  be  ascribed  to  changes  in  R(t)  [Eqs.  (4-1)  and  (4-2)]. 

The  frequency  axis  of  the  PSD  function  can  be  scaled  in 
velocity  since  v * Af/2.  One  ambiguous  velocity  interval  spans 
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the  Nyquist  range  from  -l/2At  to  l/2At,  such  that  Vamb  ■ 

X (2/2At)/2  ■ x/2At.  An  object  may  be  seen  to  fold  over  as  the 
range  rate,  which  is  measured  modulus  moves  from  one 

ambiguous  interval  to  another. 

This  effect  is  shown  in  Fig.  4-1  where  the  range  rate  of 
an  object  observed  by  a UHF  radar  gradually  changes  in  time 
(X  « 0.69  m,  At  » 1/160  sec).  Here,  consecutive  data  sets  of 
16  pulses  are  transformed  with  a conventional  FFT.  A zero 
extension  to  2S6  samples  begins  at  sample  17  and  a rectangular 
window  is  used.  There  is  a 501  overlap  of  data  from  line  to 
line,  8 old  pulses  dropped  and  8 new  pulses  added  each  time  a 
new  PSD  is  calculated  and  plotted. 

The  main  lobe  of  each  power  spectral  density  function 
locates  the  range  rate  of  the  object  during  the  processing 
interval  T.  Here  T * 0,1  second.  The  third  dimension,  a 
linear  scale  from  0 to  1,  shows  the  sidelobe  structure  of  the 
sin(TrfT)/ (tt£ T)  window  spectrum  that  has  been  convolved  with  the 
true  spectrum  (a  delta  function  at  the  true  range  rate). 

We  can  extend  each  set  of  16  pulses  to  48  pulses  with  a 4- 
point  linear  prediction  filter  before  Fourier  transformation, 
as  if  we  were  processing  0.3  second  of  data  instead  of  0.1  sec- 
ond. Figure  4-2  shows  a marked  reduction  in  sidelobe  levels 
and  a proportional  narrowing  of  the  main  lobe.  A zero  extension 
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to  256  samples  and  a rectangular  window  are  used,  although  the 
zero  extension  here  begins  at  sample  49.  A 50%  overlap  of 
original  data  points  allows  Figs.  4-1  and  4-2  to  be  compared 
line-by-line . 

Figure  4-3  displays  the  analogous  MESA  results.  The 
sharpness  of  the  peaks  and  the  absence  of  sidelobes  are  two 
immediately  apparent  features.  As  we  have  seen  earlier,  the 
error  in  locating  the  position  of  each  peak  is  considerably 
less  with  MESA  than  with  the  conventional  techniques.  Rather 
than  computing  the  full  MESA  PSD,  the  "largest  zero"  procedure 
suggested  in  Section  3.8  could  have  been  used  to  locate  the 
peaks  in  Fig.  4-3  with  the  same  precision. 

Even  when  the  data  set  is  short  and  only  a few  radar 
pulses  are  processed,  conventional  Fourier  techniques  may 
perform  satisfactorily  if  there  is  only  one  velocity  to  be 
estimated.  However,  when  the  data  are  limited,  the  conventional 
techniques  are  less  able  to  detect  closely- spaced  velocities  of 
multiple  objects  unresolved  in  range.  For  this  application  the 
linear  prediction  and  MESA  algorithms  may  prove  useful. 

4.3  Range-Doppler  Sizing 

The  size  estimation  and  imaging  of  a hard  body  rely  on 
cross-range  (Doppler)  measurements  made  in  each  of  the  range 
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cells  that  overlay  the  target.  For  example,  if  an  object  is 
spinning,  then  each  frequency  component  f^  (that  is,  each 
velocity  component  v^)  in  the  radar  return  is  proportional  to 
the  distance  r^  from  the  spin  axis  at  which  a scattering  center 
on  the  target  is  located.  That  is, 


vi  ■ *V2  ■ “spinrisinn  (4'4) 

where  X is  the  radar  wavelength,  ^Sp^n  is  the  spin  angular 
frequency,  and  0,  is  the  aspect  angle  between  the  spin  axis  and 
the  RLOS.  We  seek  the  frequency  components  f^.  With  a knowledge 
of  wSpin  and  Q,  a range-Doppler  image  can  be  constructed  from  a 
collage  of  the  cross-range  plots  from  all  range  cells  containing 
the  object. 

Figure  4-4  displays  the  evolution  of  the  conventional 
Fourier  PSD  for  the  radar  return  coming  from  the  base  of  a 
spinning  and  precessing  conical  target.  The  wavelength  X is 
5.3  cm,  At  is  0.01  sec,  and  32  pulses  are  processed  at  once 
with  a 50%  overlap  of  pulses  from  line  to  line.  Because  the 
spin  period  is  0.5  second  and  0.32  second  of  data  is  being 
transformed  per  line,  the  velocity  spectrum  is  not  time  station- 
ary and,  therefore,  becomes  smeared  within  a band  of  frequencies 
through  which  many  scattering  centers  move  during  the  processing 
interval.  The  band  is  modulated  by  the  precessional  motion, 
which  changes  the  aspect  angle  in  time.  However,  the  base 
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I radius  can  be  estimated  with  Eq.  (4-4)  and  a knowledge  of  the 

| spin  rate  and  time-dependent  aspect  angle.  The  edges  of  the 

precessional  envelope  correspond  to  the  horizons  of  the  base 
that  are  spinning  toward  and  away  from  the  radar,  so  a lower 
limit  on  the  base  radius  can  be  made  with  Eq.  (4-4)  where  the 
velocity  spectrum  has  its  greatest  width. 

Analogous  to  Fig.  4-4  is  Fig.  4-5,  for  which  the  32  pulses 
are  extended  to  96  before  Fourier  transformation.  Because  the 
data  are  not  initially  time-stationary  owing  to  the  long  process- 
ing interval,  there  is  not  a significant  improvement.  Use  of 
the  MESA  procedures  in  Fig.  4-6,  however,  does  suppress  the 
sidelobes  and  allows  a clearer  definition  of  the  precessional 
envelope. 

A reduction  in  the  length  of  the  data  set  makes  time- 
stationarity  approachable.  When  the  processing  interval  is 
small,  the  scattering  centers  on  the  target  base  cannot  exhibit 
as  large  a change  of  range  rate  and  hence  severely  smear  the 
velocity  spectrum.  However,  conventional  Fourier  techniques 
have  poorer  resolution  as  the  data  set  becomes  shorter;  in  the 
effort  to  detect  single  frequency  compon?nts  in  approximately 
time-stationary  data,  we  paradoxically  lose  the  ability  to 
resolve  them. 

J 


For  example.  Fig.  4-7  shows  conventional  Fouriepr  spectra 
when  only  8 pulses  are  processed  at  a time,  as  opposed  to  the 
32  pulses  used  earlier  (a  50%  overlap  is  retained).  Even 
though  Fig.  4-7  is  centered  where  the  precessional  envelope 
necks  down  (about  t = 7.5  seconds  in  Fig.  4-4),  there  is  no 
distinct  indication  that  the  target  is  precessing. 

If,  however,  the  8 pulses  are  linearly  extended  with  a 3- 
point  filter  to  24  pulses  before  Fourier  transformation,  the 
precessional  envelope  becomes  evident  in  Fig.  4-8.  In  this 
example,  linear  prediction  has  provided  the  same  information 
about  the  dynamic  motion  of  the  body  with  one- fourth  the  amount 
of  data  previously  used. 

Analogous  MESA  spectra  are  shown  in  Fig.  4-9.  The  absence 
of  sidelobes  more  clearly  reveals  the  precessional  envelope. 
Moreover,  there  are  suggestions  of  the  paths  of  individual 
scattering  centers  where  the  spectral  peaks  sweep  diagonally 
from  right  to  left  as  time  increases  (see  arrows).  An  ability 
to  track  individual  scattering  centers  would  be  sufficient  to 
estimate,  for  example,  the  spin  rate  of  a target  unambiguously. 

On  the  one  hand,  any  of  the  three  methods  of  Doppler 
processing  seems  to  provide  essentially  the  same  spectral 
information  if  the  data  set  is  long  or  if  it  is  not  time- 
stationary. The  spectra  computed  with  the  linear  prediction 
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and  MESA  algorithms  are,  however,  less  cluttered  with  sidelobes 
to  allow  better  definition  of  the  qualitative  features  in  the 
spectra. 

On  the  other  hand,  the  linear  prediction  and  MESA  algorithms 
can  analyze  short  data  sets  with  more  success  than  can  conven- 
tional Fourier  techniques.  Features  which  the  latter  may  fail 
to  resolve  may  be  detected  by  the  alternative  techniques. 


4.4  Other  Radar  Applications 

Only  two  applications  have  been  mentioned  here  that, 
nevertheless,  demonstrate  the  utility  of  predictive  deconvolution 
concepts.  Other  applications  might  include 

(1)  Radar  metrics  — updating  Kalman  filters 
with  fewer  radar  pulses; 

(2)  Drag  measurements  — measuring  target 
range  rates  in  less  time  and  with  more 
precision; 

(3)  Discrimination  in  clutter  - detecting 
objects  in  velocity  space  that  scatter 
weakly  compared  to  the  background  noise; 

(4)  Pulse  compression  — improving  range 
resolution  when  frequency  data  are 
digitally  transformed  into  the  time  or 
range  pulse  shape. 

This  list  is  by  no  means  exhaustive.  The  potential  user 
must  decide  if  the  linear  prediction  and  MESA  procedures 
afford  an  advantage  for  his  own  application.  In  any  case, 
these  new  techniques  are  well  worth  trying  and  offer  an  alterna- 
tive to  the  conventional  Fourier  transform. 
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Fig. 4-4.  Doppler  history  of  base  range  cell  of  spinning 
and  precessing  cone  using  conventional  Fourier  transforma- 
tion of  32  samples. 
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Fig. 4-6.  MESA  Doppler  history  analogous  to  Fig. 4-4 
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APPENDIX  I 


CALCULATION  OF  PREDICTION-ERROR  FILTER  COEFFICIENTS 

FOR  COMPLEX  DATA 


Van  Den  Bos  (1971)  has  shown  that  MESA  is  equivalent  to  a 
least- squares  all-pole  model  of  the  data  being  analyzed.  This 
means  that  a data  point  is  predicted  by  the  weighted  average  of 
its  neighbors.  The  all-pole  model  is  also  known  as  the  autoregres- 
sive model  (Box  and  Jenkirs,  1970).  Anderson  (1974)  used  this 
formalism  to  develop  algorithms  for  calculating  the  prediction 
filter  coefficients  for  a real  input  time  series.  Since  a complex 
time  series  is  often  of  interest,  we  shall  modify  Anderson's  work 
here  and  extend  it  to  complex  form.  The  same  notation  is  used 
for  ease  in  comparison.  The  reader  should  note  that  Anderson 

defines  his  filter  coefficients  as  — a„  „ instead  of  +am  „ (the 

m,n  m,n 

latter  is  the  convention  in  most  MESA  literature).  Then,  the 
MESA  power  spectrum  is  written 


P(f)  = 


P At 
m 


(1-1) 


1_  2 am,n  e 


2iTfnAt 


The  frequency  f is  limited  to  the  Nyquist  range  |f|  £ (l/2At). 
Pm  (the  residual  error  power  remaining  after  an  (m+1) -point 


v •*"* 


filter  is  convolved  with  the  data)  and  each  filter  coefficient 
a (the  nth  coefficient  out  of  the  total  m coefficients)  are 

III  | ii 

determined  by  the  solution  of  the  equation 


where  <J>(1)  is  the  value  of  the  autocorrelation  function  at  lag  1. 
For  1=0,  we  have 


N 


£ 

t=i 


xtxt 


d-3) 


as  the  (real)  variance  of  the  complex  series  (xt)  of  N points. 

'The  solution  of  (1-2)  involves  the  stepwise  increase  of  the 

matrix  dimension  by  one  and  the  determination  of  (m+2)  unknowns 

as  the  known  autocorrelation  function  is  being  calculated.  These 

unknowns  include  the  m filter  coefficients,  the  next  value  of  the 

autocorrelation  function,  and  the  error  power:  (a  a , 

r m,l  m,m 

<h (m) , P ) . There  are  only  (m+1)  equations  in  Eq.  (1-2),  so  an 
additional  relation  is  necessary.  Burg  (1968)  suggested  that  the 
additional  relation  is  the  minimization  of  the  total  error  power 
(the  sum  of  the  forward  error  power  and  backward  error  power) . 


prediction-error  filter  (1,  ^)  are 


1 1 “ J. 

"=TT  ^ lxt+i  ~ aifixtl 


d-4) 


! * 

— ry  E l*t  “ ai,ixt+i 


(1-5) 


? 2 

Minimizing  the  sum  e£‘  + with  respect  to  a^  ^ gives 


1 , 1 = 2 S Vfl/  £ Cxtxt  + xt  + lxt  + l)- 

t-1  / t«l 


(1-6) 


In  general,  for  an  (m+1) -point  filter  (1  with  the  m coefficients), 


1 

CN  — m] 


N-m 

E x! 

t=l 


am,kxt+m-k 

k=l 


(1-7) 


, N-m  m A 

(N  - m)  E xt  ~ E am,kxt+l 


d-8) 


2 2 

The  sum  e£  + e^  is  minimized  with  respect  to  am  Adding  a 
new  coefficient  will  require  editing  the  old  coefficients  (k=l, 


...  , m-1)  by  the  rule 


am,k  ” am-l,k  am,mam-l,m-k 


d-9) 


to  update  the  filter.  Equation  (1-9)  is  the  result  of  the  Levin- 
son recursion  relations  (Wiener,  1949;  Robinson,  1967)  for  the 


- . 


solution  of  Toeplitz  matrix  equations  like  Eq.  (1-2).  The  updated 
error  power  becomes 


^m  ” ^m-1^  am,mam,m^' 

If  we  set  am,o  * -1  and  am,k 
Eqs.  (1-7)  and  (1-8)  as 


(1-10) 

= 0 for  k > m,  we  can  rewrite 


2 N-m 

2 "(N  - m)  £ 
t = l 


m 


a_  . x 


k=o 


m,k  t+k 


d-11) 


'b  2 (N  - m) 


N-m 

£ 

t=i 


m 


k=o 


am, kxt+m-k 


(1-12) 


Inserting  Eq.  (1-9)  into  Eq.  ( I - 11)  and  Eq.  (1-12)  we  obtain 


N-m 


f " YJW 


— m)  ^ ^m.t  am,m^m,t 


t=l 


(1-13) 


N-m 


2 i n in 

eb  2 (N  — m)  £ ^m,t 


t = l 


a b 1 „ 

m,m  m,t 


(1-14) 


where  we  have  defined  the  series 


m 


m 


m 


,t  = 2 am- 1 , kxt+k  = £ am-l,m-kxt+m-k  (I“1S) 


k=o 
m 

bi.t  ■ S am-l, 

k=o 


k=o 
m 


kxt+m 


-k  £ am-l,m-kxt+k* 
k=o 
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The  alternate  forms  of  b + and  b'  * are  obtained  by  changing  the 

m , t ni , t 

index  k to  m-k  and  reversing  the  order  of  summation  (commutativ- 

2 2 

ity) , Minimizing  the  sum  e£  + with  respect  to  am  m gives 


N-m 

a = 2 £ b*  j.b ' 

m,t  m,t 


m.m 


m,t^m,t  + 


K.tK.S- 

(1-17) 


It  is  easy  to  show  that 

* 

b , ' b * . 2 i b « . 

m,t  m-l,t  m-l,m-l  m-l,t 


(1-18) 


^m,t  ~ ^m-l,t+l  am- 1 ,m- l^m- 1 , t+1 


(1-19) 

can  be  constructed  from  their 
previous  values  as  m increases  by  one.  The  initial  values  are 


so  that  the  series  b . and  b'  . 

m, t m, t 


and 


o,t 


= b' 


o,t 


= x. 


(1-20) 


l,t 


_ x 


t + 1 


(1-21) 


Then  b and  b'  . can  be  calculated  from  Eqs.  (1-18)  and  (1-19) 

111  ^ 111  ^ L 

as  m is  incremented  from  1 to  the  desired  filter  length  M. 

Clearly,  the  solution  of  (1-2)  is  a "bootstrap"  process 
based  on  a recursion  from  m **  1 to  M.  Figure  1-1  shows  a flow 
chart  for  the  recursive  procedure  that  calculates  the  complex 
filter  coefficients.  Following  Anderson's  Figure  1,  bm  t is  bl 


54 


M 


BWrewfflt  f- ■ 
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fr  * 


and  b'  t is  b2.  The  array  aa(t)  is  temporary  storage  of  the 
filter  coefficients  before  being  updated  into  array  a(t)  by 
Eq.  (1-9),  and  array  P(m)  is  the  error  power  updated  by  Eq.  (I- 
10)  as  m ranges  from  one  to  its  maximum  specified  value.  Element 
P(M)  is  the  final  error  power  used  in  the  MESA  spectrum. 

Figure  1-2  gives  a FORTRAN  listing  for  the  flow  chart  of 
Fig.  1-1.  The  arguments  of  the  subroutine  are 


number  of  data  samples 

complex  array  of  data  samples 

number  of  filter  coefficients  to  be 
calculated 

complex  array  of  coefficients 
real  array  of  updated  error  power 
initial  variance 

work  arrays. 


NPTS 
X 

NC0EFF 

A 
PM 
P0 
AA 

’ B1 
B2 

The  storage  for  all  arrays  must  be  supplied  by  the  calling  pro- 
gram. Inputs  are  NPTS,  X,  and  NC0EFF. 

This  subroutine  can  be  used  with  the  program  CPSPEC  given  in 
Fig,  1-3  to  calculate  a MESA  power  spectrum.  The  argument  names 
are 

NC0EFF  number  of  filter  coefficients 

A complex  array  of  coefficients 
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c 

'L 

i. 

s; 

£• 


r 

;■ 

»s 

S'* 

t. 

i' 


PSPEC  array  of  valves  of  power  spectrum 

FREQ  array  of  frequencies  at  which  PSPEC  is 

calculated 

NPSPEC  number  of  frequencies  in  array  FREQ 

PLAST  PM(NC0EFF)  from  above  = Pm  in  Eq.  (1-1) 

DT  sample  spacing  of  data  points. 

Inputs  are  NC0EFF,  A,  FREQ,  NPSPEC,  PLAST,  and  DT.  Sub- 
routine CPSPEC  is  written  with  Anderson's  negative  coefficient 
convention  (Eq.  (1-1)).  Again,  array  storage  must  be  supplied  by 
the  calling  program. 


Figure  1-4  lists  the  subroutine  LNPRED  by  which  complex  data 
may  be  linearly  extended.  All  arguments  are  input,  but  array  X 
must  be  large  enough  to  accommodate  the  extension.  Anderson's 
negative  coefficient  convention  is  incorporated  into  the  code. 

N1  original  number  of  data  points 

N2  number  of  data  points  after  extension 

X complex  data  array 

NC0EFF  number  of  coefficients 

A complex  array  of  filter  coefficients. 


i 


j 

I 


STOP 


Fig .1-1.  Flow  chart 
prediction  filter  coe 
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SUBROUTINE  COER F (NPTS,  X , NCOEr P,  E,  PH,  PO , EE,  Bl , B2) 

THIS  SUBROUTINE  CALCULATES  THE  COREL  EX  PILTER  COEFFICIENTS. 

THE  ALGOR ITHHS  USED  HERE  EES  E HODIPICETION  OP  TRE 
ALGORITHMS  DESCRIBED  BI  ENDERSON  (GEOPHYSICS,  VOL.  39, FEB.  1979) 
FOR  THE  CESE  OF  E COHPLEX  SERIES.  REEL  DETE  CEN  BE  PROCESSED  BY 
SETTING  THE  IHEGINERY  PERT  OF  THE  DETE  STORED  IN  I TO  SERO. 

INPUTS  ERE  NETS, X,NCO£FF 


NPTS  * THE  HUH  BEE  OF  DETE  POINTS  IN  THE  DETE  SET 
X > E COMPLEX  ERREY  CONTE1NING  THE  ORIGINEL  DETE 

NCOEFF  • THE  NUMBER  OF  FILTER  COEFFICIENTS  TO  BE  CELCULETED 
E • THE  ERREY  CQNTEI HI NG  THE  COMPLEX  FILTER  COEFFICIENTS 

P.1  ■ n EEL  ERREY  CONTEINING  THE  UPDETED  ERROR  PONER 

PO  « REEL  VERIEHCE  OF  THE  DETE 
EE.  Bl,  B2  ERE  MORE  ARRAYS 

STORAGE  FOR  THE  ERREYS  RUST  HE  SUPPLIED  BY  THE  CELLIHG 

PROGRAM:  X(NPIS)  , E (NCOEFF)  , EE  (NCOEFF)  ,B1(NPTS)  , B2  (NPTS)  , PH  (NCOEFF) 

ERL  THE  MINIMUM  STORAGE  REQUIREMENTS  FOR  THIS  SUBROUTINE. 


PROGRAMMED  BY  S.B.  BOWLING,  HIT- LINCOLN  LABORATORY,  FEB.  1977. 


10 


20 


21 


22 

25 


30 


35 

30 


90 

50 


COMPLEX  X.E,EE,B1,B2,XNOR,DSN,TNO 
DIMENSION  X(1),R(1),ER(1),B1(1),B2(1),PR(1) 
1N0«CHPLX(2. 0,0.0) 

PC.-0.0 

DO  10  IT»1,NPIS 

DUMMY  * X(IT)  «CQNJG  (X  (IT)  ) 

PO*  PO+  DUMMY 
PO*  PO/FLOET  (NPTS) 

NH1*NPIS-1 
El  ( 1)  *X  ( 1) 

B2 (NMl ) * X (NPTS) 

CO  20  IT-2.NM1 
Bl  (IT)  * X (IT) 

IIMl* IT' 1 
B2  (ITM 1)  *X  (IT) 

DO  SO  M*1, NCOEFF 
Ml 1*H” 1 


NMM-NPIS-M 

IF  (M  .IQ.  1)  GO  TO  25 
DO  21  IT*1  , MM1 
EE ( IT) «E(II) 

DO  22  IT-1, NMH 

Bl  (IT)  * Bl  (IT)  -CONJG  (EE  (MM  1)  ) «B2  (IT) 

B2  (IT)  - B2  (IT*1)-EE<HH1)  *81  (IT»1) 

XNOM'CMPLX  (0,0,0.0) 

DEN*CMPLX(0. 0,0.0) 

DO  30  IT«1 «NNH 

XNOB-XNON  * B2(II)*C0NJG(B1  (IT)) 

DSN-DEN  ♦ Bl (IT) *C0NJG (Bl (XT) ) * 82  (IT)  *CONJO  (B2  (IT)  ) 


IF  ( REEL(DBB)  . EJ , 0.0)  GO  TO  35 
B (8)  * TNO*(XBOH/DBB) 

GO  TO  36 

E(M) -C1PLX(0. 0,0,0) 

PONEH-PO 

IKS  .AT.  1)  PONSR»PM  (M-1) 

DUMMY* A (H)  -CONJG  (A  (H)  ) 

PN(N)  -POWER*  (1.0  - DUMMY) 

IF  (M  .sg.  1)  GO  TC  50 
DO  90  IT-1 ,HH1 

E ( IT)  * EE  (IT)  'E  (M)  -CONJG  (EE  (H'XT)  ) 

CONTINUE 

1ETURN 

END 


1 RHP-122  (1-271 


i 


Fig. 1-2. 


FORTRAN  listing  for  flow  chart  in  Fig.I-1. 
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SUBROUTINE  CPSPEC  (NCOSFP,  A , PS  PBC,  FREQ,  NPSPIC, HAST, DT) 

THIS  SUNMCUTI ME  COMPUTES  THE  POM Ek  SPECTRUM  PSPRC  AT  THE  FREQUENCIES 
STORED  IK  FR*Q.  IT  USES  RECURSION  REUTIOMS  FOR  COS  (N*THtt  A) 

AND  SIM (M'THETR)  TO  SAVE  TINS.  IT  IS  CAST  IM  CONPLEE  FORM 
SINCE  THE  PREDICTION  COEFFICIENTS  MRS  GENERALLY  COHPLBI. 

INPUTS  ARE  NCOEFF,A,rtiig,NFSPRC,PLAST,DT. 

NCOErr  • NUMBER  OF  FILTER  COEFFICIENTS 
A » COMPLEX  ARRAY  OF  riLTER  COEFFICIENTS 

FREQ  • ARRAY  OF  FRlgU  ENCIES  AT  HHICH  THE  PONER  SP3CTBUR 

IS  TO  &E  CALCULATED 

PSPEC  • ARRAY  OF  VALUES  OF  TUE  PONER  SPECTRUM 
N PS  PEC  * NUMBAR  OF  FREQUENCIES  AT  NHICH  THE  PONER  SPECTRUM 
IS  TO  BE  CALCULATED  • DIMENSION  OF  FREQ  AND  PSPEC 
PL  AST  ■ RESIDUAL  ERROR  PONER  AFTER  THE  FILTER  HAS  OPERATED 

ON  THE  DATA  » SLEHSHT  PN(HCOSFF)  FROM  SUBROUTINE  •COEFF* 
DT  • SPACING  OF  THE  DATA 

STORAGE  FOR  ARRAYS  IS  SUPPLIED  RY  THE  CALLXRG  PROGRAM.  MINIMUM 
STORAGR  REQUIREMENTS  ARE  A(NCOIFP)  , PSPEC  (NPSP  AC)  , FREQ  (MPSPRC) 

FOR  THIS  SUBROUTINE. 

PROGRAMMED  BY  S.B.  30NL1HG,  HIT -LINCOLN  LABORATORY,  FEB.  1977. 

COMPLEX  A , ZK, DE  N 

DIMENSION  A(1)  .PSPEC  (1)  ,PREg(1> 

FACTOR*  -6.2831053*01 
Do  100  1*1 , NPSPEC 
CM- 1.0 
SN*G.O 

CEN«CHPLX(1. 0,0.0)  — 

ARG*FACTCR*FREQ<I)  _ . »■  A Hi  L 

xvSB*.  oca 

UHP*  C*HG»CN  - SAFG*SN  lj  J 1 * * 

SN-  CAFG’SN  ♦ SARG*C N 

CN»  Tssr 

IN*  CMrtX(CN.SN) 

DEN-  DIN  - A { N)  *IN 
50  CONTINUE 

ESQ*  DFN*CONJG<Oi.N| 

IF  { DSQ  . EQ.  0.0  l DSQ*1.0E-11 
PSPEC(I)*  (1.0/OSg)*PLAST*DI 
00  CONTINUE 

RETURN  

END  I RHP-1 22  (1-3) 


COM 


Fig. 1-3.  FORTRAN  listing  to  calculate  MESA 
power  spectrum. 
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SUBROUTINE  LNPMD(N1,N2,I,NC05FF,4> 

c 

C THIS  SUBROUTINE  LINEARLY  EXTENDS  Till  CONPLEX  DATA 
C FROM  N1  POINTS  TO  M2  POINTS.  THE  ORIGINAL  DSTN  , POSXTIONND 
C IN  THZ  FIRST  N 1 ELEMENTS  OP  ARMEY  X,  ERE  SHIFTED  TO  THE  MIDDLE  OF  X. 
C FOR  NERD  AND  BACKYARD  PREDICTIONS  ERE  DONE  UNTIL  THE  TOTAL 
C NUN  BAR  OF  CATE  POINTS  IS  N2. 

C 

C NOTE  THAI  THE  N2  POINTS  CONTAIN  THE  ORIOINEL  Ml  POINTS. 

C 

C INPUTS  EPS  N1 ,N2.X,NCOEFF,R. 

C 

C ARRAY  X IS  MODIFIED  ; THE  ORIOINEL  DATE  POINTS  HATE  BEEN 
C SliIF.ED  TO  THE  MIDDLE  END  PB5DICTI0NS  ERE  DONE  ON  BOTH  SIDES. 

C 

C THE  FILTER  COEFFICIENTS  SHOULD  ALREADY  HEYE  BEEN  CALCULATED 
U BY  PROGRAM  • CO I Ft. • 

C 

C N1  > ORIGINAL  NUNPRR  OF  POINTS  IN  EEREY  X 

C N2  « NUMBER  OF  PC1NTS  TO  MHIGtl  ARRAY  X IS  TO  BE  EXTENDED 

C X • COMPLEX  ARRAY  OF  DATA  SAMPLES 

C NCOEFF  ■ NUMBER  OF  PREDICTION  FILTER  COEFFICIENTS 

C A • COMPLEX  AhRAY  OF  FILTER  COEFFICIENTS 

C 

C STORAGE  FOR  THE  ARRAYS  SHOULD  BE  SUPPLIED  BY  THE  CALLING 
C PROGRAM:  MINIMUM  STORAGE  RFg'J  IRSSENTS  ARE  X (N21 , A (NCOEFF) 

C FOR  THIS  SUBROUTINE. 

C 

C PROGRAMMED  BY  S.  B.  ROLLING.  NIT-LINCOLN  LABORATORY,  FEB.  1977. 

C 


CjHPLEX  l.k 
DIMENSION  X(1),A(1) 


C 

c 


c 

c 

c 


SET  UP  LIMITS  FOR  DO  LOOPS 

L 1*  N 2/2  - N 1/2 

L2«  B2/2  ♦ N1/2 

IF | MOD(Nl.i)  .EO.  1 I L2*L2»1 

SHIFT  ORIGINAL  PAT  A TO  MlDDlt  OF  ARRAY  X 


CO  IOC  1*1, N1 
J=  NT  - (1-11 
K*  L2  - C-1) 
100  X(K)*X(J) 


C DO  FOKRARL  PREDICTION 
C 

NS*  N2-L2 
CO  *00  I* 1 , N 3 
J*  L2*I 

X (0)  *CMPLX  (0. 0, 0. 0) 

DO  20C  K-1. NCOEFF 
ECO  X(J|*  X(J)  ♦ A(K)«X(J-K) 

C 

C DO  BhCKHARC  PPLDICTION 

c 

EO  SOC  1*1, LI 
J»  LI-  (1-1) 

X(J)*CMPLX  (w.O.C.O) 

CO  SOO  K*1  # NCOEFF 

J00  X(J)»X(J)  + CONJG  (A(K)1*X(J*K) 

RETURN  I RHP-1 22  (I-4)~ 

IND 


c 


Fig. 1-4.  FORTRAN  listing  to  linearly  extend  a complex 
data  set. 
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APPENDIX  II 


A COMPARISON  OF  SOME  OF  THE  CHARACTERISTICS 
OF  MESA  AND  CONVENTIONAL  FOURIER  SPECTRAL  ANALYSIS 

Here  we  tabulate  some  of  the  salient  features  of  the  two 
methods  of  spectral  analysis.  For  a more  mathematical  presen- 
tation, the  reader  is  referred  to  Lacoss  (1971)  and  to  Chen  and 
Stegun  (1974).  The  notation  in  the  table  is 

N = number  of  values  of  known  autocorrelation 

function 

P * power  of  a pure  tone  in  the  spectrum 

At  * spacing  of  input  data 
NPTS  a number  of  data  samples. 
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COMPARISON  OF  MESA  AND  FFT  CHARACTERISTICS 


APPENDIX  III 


ON  INTEGRATING  THE  MESA  POWER  SPECTRUM 
TO  FIND  THE  POWER  SPECTRAL  DENSITY 

Lacoss  (1971)  has  pointed  out  that  it  is  the  MESA  power 
spectral  density  which  is  the  more  appropriate  measure  of  the 
relative  power  of  spectral  components.  Since  the  power  spectrum 
P(f)  is  calculable  at  any  frequency  |f|  £ l/2At  (assuming  the 
data  are  band-limited  and  properly  sampled),  in  principle  it  is 
easy  to  perform  a numerical  integration  over  a small  bandwidth  6f 
to  calculate  a power  spectral  density.  However,  some  of  the 
peaks  in  P(f)  may  be  so  narrow  and  sharp  that  they  may  be  diffi- 
cult to  detect  if  the  initial  grid  spacing  in  frequency  is  coarse. 

Radoski  __t  al . (1975)  emphasize  that  if  the  signal- to-noise 
ratio  is  high,  a coarse  frequency  grid  can  be  insufficient  to 
determine  the  actual  locations  of  spectral  peaks  in  the  MESA 
power  spectrum.  As  the  number  of  filter  coefficients  increases, 
the  peaks  tend  to  approach  line  spectra  (delta  functions). 

Radoski  trt  al . (1975)  have  suggested  a systematic  search  procedure 
(used  in  this  report  in  a modified  form)  to  calculate  the  power 
spectral  density  from  the  MESA  power  spectrum: 

(1)  The  data  are  normalized  to  unit  variance.  Each  datum  x 
is  replaced  by  (x  - x)/a  where  x and  a are  the  (com- 
plex) average  value  and  (real)  standard  deviation  of 
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the  original  data.  By  Parseval's  Theorem,  the  integral 
of  the  power  spectral  density  (PSD)  must  be  unity. 
Therefore  we  can  feel  confident  that  the  PSD  has  been 
correctly  estimated,  the  spectrum  integrated,  and  all 
spectral  lines  detected,  if  numerically  J*PSD(f )df  « 1 . 

(2)  Triplets  of  points  in  the  power  spectrum  are  examined 
such  that,  for  (f^  + 1 - f^_^)  = ^f,  whenever 

PCf^P  < P(f.)  > P(fi+1)  (i.n-i), 

the  midpoint  is  near  a possible  spectral  peak.  P(f)  at 
two  intermediate  points  is  calculated  and  the  five  are 
examined  to  extract  a new  triplet.  Simpson’s  rule  is 
used  to  integrate  P(f)  as  the  search  procedure  works 
its  way  up  into  the  peak  until  two  calculations  of  the 
integral  over  the  bandwidth  6f  are  within  1%  of  each 
other.  Then  the  procedure  shifts  to  the  next  triplet 
of  points  in  the  power  spectrum  spanning  6f  and  the 
integration  process  is  repeated.  If  Eq.  (III-l)  is  not 
satisfied,  Simpson's  rule  integrates  P(f)  over  Sf  with 
enough  points  so  that  two  successive  integrations  are 
within  1%  of  each  other.  See  Radoski  e^t  al . (1975)  for 
details . 

Normalizing  the  data  to  unit  variance  removes  the  dc  or  zero 
frequency  part  of  the  spectrum,  which  may  not  be  appropriate  for 
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some  applications.  If  the  data  are  not  normalized,  then  J*PSD(f)df 

should  equal  the  variance  a of  the  data  (which  is  equal  to  the 

zero-lag  value  of  the  autocorrelation  function).  Depending  on 

the  form  of  Parseval's  Theorem  used,  JpSD(f)df  and  a2  may  differ 

1/2 

by  a factor  of  2 it  or  (2ir)  ' . (Parseval's  Theorem  states  that 
the  total  average  power  (or  mean- square  value)  of  x(t)  is  equal 
to  the  integral  of  the  power  spectral  density  over  all  relevant 
frequencies. ) 
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